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Abstract 


In  an  earlier  paper,  approximate  equations  of  motion  were  derived 
which  are  applicable  to  nondissipative , very  nonadiabatic,  buoyant  flows 
of  a perfect  gas.  In  the  present  paper,  these  approximate  equations  are 
recast  in  a form  in  which  they  can  be  integrated  by  finite  difference 
techniques.  These  nonlinear  equations  are  made  dimensionless,  specialized 
to  two  dimensions  and  linearized.  Finite  difference  approximations,  using 
central  differences  in  space  and  leapfrog  in  time,  are  made  to  the  linear- 
ized equations.  To  start  the  computation,  a first-order  scheme  is  used 
for  the  initial  time  step.  The  dependent  variables  are  density,  the 
horizontal  and  vertical  velocity  components  and  pressure,  and  these 
variables  are  defined  at  various  positions  within  a grid  cell  (see  Figure 
lb) . The  computer  time  required  for  various  size  computations  and  the 
accuracies  obtained  are  reported. 

At  each  time  step,  a large  sparse  system  of  linear  algebraic  equations, 
the  finite  difference  approximation  to  the  elliptic  equation  for  the  pres- 
sure, must  be  solved.  A hybrid  method,  combining  the  interative  algorithm 
called  conjugate  gradients  with  a fast  direct  method  for  solving  Poisson's 
equation,  is  used  to  solve  the  algebraic  system  efficiently  and  accurately. 
(Typically,  the  solution  of  the  algebraic  equations  arising  on  a 31  x 31 
mesh  is  determined  to  five  significant  figures  after  between  two  and  five 
iterations  of  the  conjugate  gradients  algorithm  using  from  two  to  five 
seconds  of  CPU  time  on  the  N.B.S.  UNIVAC  1108.) 

Finally,  results  of  some  computations  are  presented  using  computer- 
generated plots.  A discussion  is  given  of  the  software  required  to  generate 
these  plots.  A very  brief  discussion  is  presented  of  the  analytical  pro- 
cedures required  to  obtain  a solution  to  the  differential  equations  describing 
heating  of  a homogeneous  fluid.  Computations  of  the  solution  to  the  dif- 
ference equations  are  found  to  confirm  the  conclusions  of  the  analysis. 

Also,  these  computations  are  found  to  yield  flow  fields  which  behave 
qualitatively  the  same  as  flow  fields  observed  in  a room  fire.  Examples  of 
computations  of  the  flow  fields  induced  by  heating  in  a stratified  fluid 
are  also  presented  and  discussed. 
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I . Introduction 

Convection  in  fluids  heated  from  below  has  been  of  interest  for  over 

one  hundred  years.  Much  of  this  interest  has  arisen  from  astrophysical 
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and  geophysical  research.  More  recently,  technological  applications  in 
areas  such  as  nuclear  reactor  design,  liquid  natural  gas  transport,  and 
fire  safety  have  led  to  the  analysis  of  phenomena  occurring  injvery  dif- 
ferent ranges  of  geometrical  and  physical  parameters  than  those  considered 
in  earlier  studies. 

The  principal  application  of  interest  to  the  authors  is  buoyant 
convection  induced  in  room  fires.  The  National  Bureau  of  Standards, 
acting  through  the  Center  for  Fire  Research,  is  involved  in  a wide  rante 
of  investigations  into  various  pjiysicaland  chemical  phenomena  arising  in 
unwanted  fires.  The  buoyant  convection  studies  reported  here  are  part  of 
that  effort.  The  heavy  mathematical  and  computational  emphasis  in  this 
project,  together  with  its  intrinsic  interest,  have  led  to  its  sponsor- 
ship as  a joint  activity  with  the  Center  for  Applied  Mathematics.  This 
report  is  meant  to  record  the  current  status  of  the  project.  As  such,  it 
should  be  considered  an  interim  working  document,  rather  than  an  account 
of  a completed  piece  of  work. 

The  approximate  equations  governing  the  large  scale  motion  of  a gas 

in  an  enclosure  driven  by  a prescribed  volumetric  heat  source  were  derived 

4 

in  an  earlier  paper.  The  resulting  mathematical  model,  a ’’thermally  ex- 
pandable fluid,"  is  rather  different  from  other  idealized  fluid  equations, 
such  as  the  Boussinesq  model,  usually  employed  in  buoyant  convection 
studies.  The  differences  arise  from  the  need  to  consider  large  density 
variations  in  a closed  geometry  of  limited  vertical  extent.  This 
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combination  of  parameters  also  arises  in  some  boiling  water  nuclear  reactor 

studies,  where  a single  fluid  approach  to  the  two  phase  flow  results  in  a 

22  21 

thermally  expandable  fluid  model.  * 

Due  to  the  complexity  of  the  equations  of  interest  and  the  limited 
experience  in  dealing  with  them,  considerable  care  was  necessary  in  de- 
veloping numerical  computation  techniques.  The  thermally  expanding  fluid 
motions  share  some  common  features  with  both  incompressible  fluid  flow 
and  with  the  hydrodynamic  models  used  in  numerical  weather  prediction. 

A sampling  of  the  techniques  used  in  the  analysis  of  these  related 
equations  may  be  found  in  references  5-9.  A study  of  these  methods  led 
the  authors  to  the  conviction  that  the  computational  procedures  should 
be  developed  on  a related  system  of  linear  equations  before  nonlinear 
calculations  were  attempted. 

Three  features  of  the  thermally  expanding  fluid  model  which  were 
felt  to  be  of  particular  importance  are  amenable  to  analysis  within  the 
context  of  a linear  problem.  They  are  the  numerical  computation  of 
internal  wave  oscillations  over  many  cycles,  efficient  calculation  of 
the  self  adjoint  elliptic  equation  governing  the  pressure  distribution, 
and  the  use  of  a distributed,  volumetric  heat  source  of  prescribed  strength 
(as  opposed  to  boundary  conditions)  to  drive  the  motion  in  a rectangular 
enclosure.  The  linear  system  was  chosen  by  formally  linearizing  the  thermally 
expanding  fluid  model  about  a rest  state  stably  stratified  in  the  vertical 
direction.  The  linearized  motion  is  driven  by  the  prescribed  heat  source. 

The  homogeneous  solutions  to  the  problem  posed  above  were  studied  by 
the  authors  in  reference  10.  These  solutions  are  internal  waves  which 
arise  due  to  the  stable  stratification.  These  waves  were  used  to  test 
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finite  difference  schemes  for  accuracy  and  stability.  Second  order  accurate 
spatial  differences  were  used  in  conjunction  with  both  first  and  second 
order  accurate  time  differencing  to  derive  candidate  finite  difference  approx- 
imations to  the  homogeneous  (internal  wave)  equations.  The  finite 
difference  equations  were  solved  exactly  by  analytical  means  and  compared 
with  the  known  solutions  to  the  continuous  problem.  The  finite  difference 
scheme  used  in  the  present  study  to  solve  the  inhomogeneous  equations  was 

selected  on  the  basis  of  the  analysis  reported  in  Reference  10. 

This  report  describes  the  numerical  implementation  of  the  finite 

difference  scheme  selected  on  the  basis  of  the  analysis  outlined  above. 

The  motions  are  driven  by  a prescribed  distributed  heat  source  located  at 
the  bottom  of  a rectangular  enclosure.  Section  II  presents  the  non-linear 
continuous  equations  derived  in  reference  4 and  describes  the  linearization. 
The  finite  difference  equations,  incorporating  second  order  central  dif- 
ferences in  space  and  leapfrog  in  time  are  shown  in  section  III.  The 

solution  of  the  pressure  equation,  based  on  fast  direct  methods  for  solving 
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Poissons  equation  combined  with  an  iterative  conjugate  gradient  technique  ’ 

is  discussed  in  section  IV.  Section  V outlines  the  solution  procedure  for 
the  coupled  set  of  linear  finite  difference  equations.  Several  examples 
of  computational  results,  together  with  the  computer  graphics  used  to 
display  them,  are  described  in  Section  VI.  Finally,  a summary  of  the  work 
presented  in  ■ this  report  is  given  in  Section  VII. 
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II.  FORMULATION  OF  CONTINUOUS  PROBLEM 

A.  Full  Nonlinear  Equations 

4 

In  an  earlier  paper  nonlinear  equations  describing  three 
dimensional  thermally-driven  buoyant  flows  in  an  enclosure  had  been 
derived  and  their  properties  discussed.  In  this  section  these 
nonlinear  equations  will  be  rewritten  in  a form  appropriate  for 
numerical  integration  by  finite  difference  techniques.  This  recast 
set  of  equations  is  then  linearized  and  the  boundary  conditions  for 
the  linear  equations  are  presented. 

As  in  Reference  4 we  consider  an  inviscid,  non-heat-conducting 
perfect  gas.  The  magnitude  and  the  spatial  variation  of  the  heat 
source  (representing  the  exothermic  reaction  in  a fire)  are  taken 
as  known;  justification  for  such  a model  is  given  in  Reference  4 . 

The  fluid  and  the  fire  source  are  assumed  confined  in  a closed 
rectangular  room  with  the  center  of  the  source  along  the  floor.  In 
contrast  to  Reference  4,  we  consider  only  a completely  enclosed  room  (no 
leaks),  and  we  confine  attention  in  the  linear  problem  to  the  two  di- 
mensional evolution  of  the  flow. 

Equations  (ll)  of  Reference  4 are 


fo  ft)  = p£T 


( i ) 


Here 


is  density,  the  velocity  in  the  i — coordinate  direction 
(i  » 1,  2,  3),  p is  the  pressure  excess  above  the  mean  pressure 
PQ(t)  in  the  room,  T the  temperature,  the  constant  pressure 
specific  heat,  R the  gas  constant,  is  the  gravitational 
acceleration  and  Q(x^,  t)  the  specified  volumetric  heat  source. 

The  spatially  uniform  mean  pressure  pQ(t)  depends  only  upon  time 
and  increases  because  of  the  heating  within  the  room.  It  is 
determined  in  a completely  enclosed  room  by  the  equation 

^ = Vjv  <?(*••>■«  4V  M 

where  ^is  the  ratio  of  specific  heats,  V is  the  volume  of  the  room 
and  the  integration  is  performed  over  this  entire  volume.  Equation 
(2)  is  a thermodynamic  statement  that  the  mean  pressure  rise  as  a 
function  of  time  is  determined  by  the  total  heat  added  to  the  room. 

(Heat  can  only  be  added  or  removed  volumetrically  and  not  through 
the  walls  because  thermal  conduction  and  radiative  transport  have 
been  ignored  in  this  model.) 

Equations  (l),  derived  in  Reference  4,  are  an  approximate  set 
of  nonlinear  equations  which  are  applicable  to  highly  nonadiabatic, 
nondissipative  buoyant  flows  of  a perfect  gas  in  three  dimensions. 

These  equations  were  derived  under  the  assumption  that  heat  is  added 
slowly  compared  with  the  time  required  to  equilibrate  the  pressure 
over  the  spatial  extent  of  the  heat  source,  an  assumption  appropriate 
to  flows  induced  by  fires.  They  are  characterized  by  the  fact  that 
the  spatially  uniform  mean  pressure  appears  in  the  energy  and  state 
equations  while  the  spatially  variable  portion  of  the  pressure  appears 
only  in  the  momentum  equation.  The  equations  admit  buoyant  or  internal-wave 


6 


motions  while  "filtering  out"  high-frequency,  acoustic  waves.  They 
reduce  to  the  Boussinesq  equations  when  heating  is  mild,  total  density 


variations  are  small,  and  variations  in  the  mean  background  pressure  can 
be  neglected  (as  would  be  the  case  if  the  room  considerated  here  were  open 
or  if  the  mean  pressure  variation  were  comparable  to  the  spatial  pressure 


perturbation.) 

Equations  (l),  fully  nonlinear  and  three-dimensional,  will  be 
recast  into  a form  suitable  for  numerical  computation.  To  do  this, 
we  take  the  substantial  derivative  of  the  equation  of  state  and  use 
this  with  the  energy  equation  to  eliminate  the  temperature.  The 
resulting  equation  describes  the  evolution  of  the  density  under  heating 


where 


= Tpo  - iijr] 


(3) 


(+) 


Equation  (3)  and  the  continuity  equation  identify  D(x^ » t)  as  the  divergence 

2Ui  _ 

Finally,  as  in  Reference  4,  the  equation  for  the  spatially  variable 


D (xz , ir) 


(sj 


portion  of  the  pressure  is  obtained  by  dividing  the  momentum  equations 
by  density  and  taking  the  divergence  of  these  equations.  The  resulting 
equation  is  r~ 

^k) =~[^r(a 

The  boundary  conditions  on  these  equations  are  that  velocity  normal 
to  any  (impermeable)  wall  vanish. 

Hi  /ni  = o 


'dx-i  V P 


^3  U-c 


f) 

A / 


Mr 


u 


rn 


where  are  the  normal  components  of  a vector  describing  the  boundary 
walls.  From  Eqs.  (l)  and  these  conditions, the  appropriate  boundary 
conditions  on  the  pressure  equation  are  obtained 
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*‘Tfc-  “ P3'n:*i 


The  complete  set  of  recast  nonlinear  equations  are  gathered  and  re- 
written below 

%f,  =-p 


(.  _ -J — S»t> 


--  !3xrK§^)+ 


■4e  ~ -V1-  Q(k^V  11/ 

= Y^fcf  QOr-D  <?(*:,■<?- 

with  boundary  conditions 

uc  /n-i  = o , ^ 

B.  Linear  Equations 

In  this  subsection  the  equations  derived  above  are  linearized 
and  specialized  to  flow  in  two  spatial  dimensions  (x  = x^,  y = Xg) 
in  a rectangular  enclosure,  0 £ x <_  L,  0 £ y <_  H.  The  fluid  in  the 
enclosure  initially  is  assumed  to  be  arbitrarily  stably  stratified 
with  density  ^(y). 

For  convenience  we  will  rewrite  variables  in  terms  of  dimensionless 
quantities.  This  procedure  can  be  applied  to  the  three-dimensional, 
fully  nonlinear  equations  also,  but  we  will  limit  ourselves  here  to 
presenting  only  the  equations  in  two  dimensions  in  linearized  form. 

Let  u be  the  velocity  in  the  horizontal  or  x-direction  and  v the 
velocity  in  the  vertical  or  y-direction.  Define  the  dimensionless 
variables,  denoted  by  hats,  as  follows: 


(*) 


0°) 
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u.=  si^T1  f-  0.  v=  e Lr 

<J>  = C (^  -+•  e.  § Q 

p-o  fo  (t)  -+  ^ £P0  (|)  + & f J (II  ) 


where 

and 


* = H £ , 

t~  (H/f)'^  t 

<?-  Q0$C**?J) 

f©  ft)  = I + £ £,  ft) 

£ = (H/g.)*.  <?.(dr-ij 


Here  ^,is  the  ambient  density  and  the  pressure  at  a specified 
altitude,  here  taken  to  he  y=o.  Qq  is  the  magnitude  of  the  volumetric 
heat  source,  and  e is  a parameter  relating  the  magnitude  of  the  heat 
source  to  ambient  pressure  (or  the  energy  per  unit  volume  in  the 
ambient  gas ) . 

For  a heat  source  characteristic  of  the  rate  of  heat  release 
within  a small  room  fire,  the  parameter  e is  often  small.  Provided  that 
the  time  does  not  become  too  large,  then,  so  that  the  cumulative  effects 
of  the  heat  addition  remain  unimportant,  e can  be  used  as  a small  parameter 
and  Eqs . (10)  can  be  formally  linearized  about  the  stably  stratified  rest 

state.  The  resulting  two  dimensional  equations  are 


(§)  £ 


(13) 

I 
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p i r\  !&  _ 

8(*,p)  - T - * 1 ^ I 0(£,$,t)] 


-^-  =*T«<f  i*«  q(*\§£) 

The  aspect  ratio  AR  of  the  room  is  the  ratio  of  the  height  H to  the 
length  L of  the  rectangular  room. 

The  boundary  conditions  corresponding  to  Eqs.  (13)  are  that  no 
mass  flow  from  the  enclosure: 


tl-  o 

II 

0 

o^-dl  )c  — , 

A 

ir  = o 

^=o 

cu»^X  =-  ! 

These  conditions,  Eqs.  (I*f),  can  be  transformed  into  corresponding 
boundary  conditions  on  the  elliptic  equation  for  the  pressure,  the 
fourth  of  Eqs.  (l3),  using  the  two  momentum  equations  in  Eqs.  (13) 


y ~ o & 


= )hz 


_ as,  _ _ 6 

^ ^ 3 ir  » 


= ( 


(W 


Os) 


In  the  computations  presented  later,  two  forms  for  the  ambient 
stratification  have  been  selected.  For  most  calculations  an  exponential 
stratification  have  been  used: 

£><£)  = (~%/ys) 

This  form  includes  the  special  case  when  Yg  is  very  large  so  that  is 

practically  constant.  For  some  of  the  computations  a "two  layer"  ambient 
stratification  was  used:  ^ 

9o(§l-«f  C-3/%) 

P.<3>  = >/ys„-  y‘/y,  ] <5  £ i 

The  heat  source  was  assumed  to  have  a simple,  separable  form: 


<y*o 


- 10  - 

where  f (t ) = tanK  a£,  A is  a constant  determining  the  rate  at  which 

A 

the  heat  source  is  "turned  on,"  and  xc  is  the  location  along  the 
"floor"  (y=o)  at  which  the  center  of  the  heat  source  is  located. 

A 

This  form  for  the  heat  source  Q was  chosen  as  representative  of  a 
fire  for  analytical  convenience  only;  another  form  could  equally  have 
"been  chosen. 

Several  important  observations  must  be  made  about  Eqs.  (l3)-(l6). 

First,  we  note  that  D(x,  y,  t),  from  the  form  given  in  Eqs.  (l3), 
integrated  over  the  rectangular  region  is  zero.  This  fact,  together 
with  the  boundary  conditions  (IS),  imply  that  the  pressure  equation 
the  third  of  Eqs.  (|3),  when  integrated  over  the  rectangular  region 
is  identically  satisfied.  Such  a condition  is  necessary  if  the  pressure 
equation  is  to  have  a solution  at  all  since  this  elliptic  equation 
has  Neumann  boundary  conditions.  When  the  pressure  equation  has  a 
solution,  this  solution  is  only  determined  to  within  a constant,  and  the 
constant  is  determined  to  be  consistent  with  the  initial  condition,  or 
constant,  chosen  in  the  integration  of  the  mean  pressure  equation,  the 
last  of  Eqs . (13) . 

Second,  it  should  be  noted  that  the  pressure  equation  in  this 
linearized  formulation  is  in  fact  separable:  multiplication  of  the 

pressure  equation  by  ^Q,(y)  shows  that  the  coefficients  of  the 
derivatives  of  y in  this  equation  are  only  functions  of  y (the  coefficients 
of  derivatives  of  x are  unity).  However,  in  the  nonlinear  formulation 
of  this  problem,  the  pressure  equation  is  nonseparable . Since  the 
linearized  equations  are  being  used  to  test  solution  techniques,  the 
solution  procedure  for  the  pressure  equation  discussed  in  Section  IV 
is  for  nonseparable  elliptic  equations. 
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III.  FINITE  DIFFERENCE  EQUATIONS 

In  this  section  the  finite  difference  equations  used  to  approximate 
the  linear  partial  differential  equations,  Eqs.  (13),  and  the  boundary 
conditions,  Eqs.  ( IH)  and  (15),  will  be  presented  and  discussed. 

These  difference  equations  are  written  on  a "staggered  mesh"  in  the  style 
of  reference  6 and  were  chosen  to  be  approximations  to  the  differential 
equations  which  are  second-order  accurate  in  both  space  and  time. 

In  Figure  la,  the  rectangular  enclosure  in  dimensionless  variables 
is  shown  together  with  a schematic  representation  of  the  spatial 
grids  used  for  the  finite  difference  scheme.  The  grid  formed  from 
solid  lines  represents  the  basic  mesh  into  which  the  enclosure  is 
divided:  in  general  there  axe  M mesh  cells  in  the  x-direction  and 

N mesh  cells  in  the  y-direction.  (The  hat  designation  on  dimensionless 
variables  will  now  be  dropped  for  convenience.  Subsequently  all  quantities 
will  be  understood  to  be  dimensionless.')  Upon  this  basic  mesh,  the  two 

components  of  the  vector  velocity  (u,  v)  and  single  component  of  the  vector 

A ^ , 2vu-_  BU 

vorticity  0>  = -jx  ^are  defined. 

The  second  grid,  formed  by  joining  the  center  points  of  the  basic 
grid  cells  and  denoted  by  dashed  lines,  is  that  upon  which  scalar 
quantities  such  as  density  (0  and  pressure  >ys>are  defined.  In  Figure  la 
the  densities  in  the  left-hand  column  of  cells  and  in  the  bottom  row 
of  cells  are  shown  to  indicate  how  they  are  enumerated  for  the  numerical 
computation. 

In  Figure  lb  a typical  mesh  cell  is  shown,  illustrating  where  all 
of  the  dependent  variables  in  the  finite  difference  scheme  are  defined 
relative  to  the  cell. 
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Fig.  lb  A typical  mesh  cell,  with  center  located 
at  x = (i  - 1/2)  6x  and  y = (j  - 1/2)  Sy,  illustrating  where  all 
dependent  variables  for  the  finite  difference  scheme  are  defined. 


lb 


The  following  discretely  evaluated  functions  will  denote  approximations 
to  the  corresponding  analytical  solutions  to  Eqs . (13): 

U.T;  = U.(lcU,  /"<*) 

d 


fij^  jo ■>  (j-fcO-hj- , 

DT;  ^ , p*,j  = 

60  ^ = 60  ( * f*.,  £ /7t  cT-6y) 

where  5 x = AR/M  and  6 y = l/N  are  the  mesh  cell  sizes  in  the  x-  and  y- 
directions  respectively  and  where  6t  is  the  time-step  size.  Such  a 
staggered  grid  is  commonly  used  for  multidimensional  finite  difference 
integrations^ 

With  this  notation,  the  following  set  of  finite  difference  equations 
was  used  to  approximate  Eqs.  (13): 
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k = mn 
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Spatial  derivatives  in  Eqs . (13)  are  approximated  by  central  differences 
and  leapfrog  is  used  for  time  derivatives.  These  equations  have  second 
order  truncation  errors  in  both  space  and  time;  i.e.,  the  local  discret- 
ization errors  made  in  approximating  the  differential  equations  (13)  by 
these  difference  equations  is  0(5x  ),  0(6y  ),  OCSt^) „ In  the  first  three 
of  Eqs.  (19),  the  time  derivatives  have  been  replaced  by  central  dif- 
ferences which  give  explicit  formulae  for  computing  dependent  variables 
at  a new  time  level,  a difference  method  commonly  referred  to  as  the 
leapfrog  method.  With  this  method  dependent  variables  at  two  previous 
time  levels  are  required  to  compute  the  quantities  at  a new  time  level; 
hence  the  scheme  is  said  to  be  a three-time-level  scheme.  A starting 
procedure  is  required  to  generate  the  dependent  variables  at  the  first 
time  level,  given  these  quantities  initially  (at  t=0) . 

To  obtain  dependent  variables  at  the  first  time  level,  an  explicit 
first-order  (in  time)  procedure  is  used: 
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where  all  supplementary  quantities  are  defined  in  Eqs.  (20). 

The  computational  procedure  by  which  Eqs.  ( 19 ) and  (2l)  are  solved 
is  very  important.  At  first  inspection,  it  may  appear  that  Eqs.  (21) 
are  implicit  for  example.  However,  when  Eg.  (21a)  is  solved  first  for 
density.  Eg.  (21d)  second  for  pressure  and  the  remaining  equations  are 
solved  subsequently,  Eqs.  (2l)  are  found  to  be  explicit.  Also,  the  time 
level  associated  with  the  discretely  evaluated  variable  representing  pres- 

oire  may  not  be  the  expected  one.  These  considerations  are  discussed  in 

Section  V. 

The  two  forms  for  the  ambient  density  stratification,  given  by  Eqs. 

(16a)  and  (l6b)  in  the  continuous  case,  have  been  used  in  discretized  form 
in  all  of  the  computations  presented  in  Section  VI.  For  most  computations 
n exponential  stratification  was  used 

C~  ( i~  &.)  ,]  (=?  * a) 

We  note  that>  when  Ys  is  taken  to  be  large  (generally  106  in  our  compu- 
tations), the  ambient  density  is  essentially  constant.  For  some  compu- 
tations, a discretized  "two-layer"  model  of  the  ambient  stratification 

was  used;  f- 

pOM-  = ^ L - Q -y*)fy/ys]  l-ji  j~s 

In  all  cases,  the  initial  velocities  Ufj  , UT ij  and  the  density  deviation 
from  ambient  were  taken  to  be  zero. 

Boundary  conditions  for  these  difference  equations  are  obtained  from 
Eqs.  (l4)  and  ( 15 ) ; they  are 

i r-  /yv-  /V1 

- ^7,  m — cJ  if^ 
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IT 


and 


$sN- 


The  pressures  and  densities  denoted  by  a tilde  are  values  defined  beyond 
jthe  grid  defined  in  Figure  1.  These  values  can  be  regarded  therefore  as 
quantities  defined  in  fictitious  cells  surrounding  the  rectangular  region 
of  interest.  However,  in  the  procedure  for  solving  the  linear  algebraic 
system  of  equations  for  the  pressure,  the  boundary  conditions  are  incor- 
porated directly  into  the  system:  the  pressure  and  density  values  in  the 

fictitious  cells  cancel  and  never  appear  in  the  final  system  of  equations 
for  the  pressure  in  the  rectangular  region. 

The  procedure  used  to  solve  the  discretized  pressure  equation  is 
discussed  in  Section  IV.  Before  this  discussion,  however,  we  note  that 
the  discrete  pressure  equation  together  with  its  boundary  condition 
retains  the  important  property  of  the  partial  differential  equation  for 
the  pressure  with  its  boundary  conditions  noted  at  the  end'  of  Section  II. 
The  definition  of  D<^  presented  in  Eqs.  (20)  shows  that 


This  fact,  plus  that  concerning  the  boundary  conditions  noted  above,  allows 
one  to  demonstrate  that,  for  the  discretized  pressure  equation,  the  left 
hand  side  and  the  right  hand  side,  when  summed  over  all  grid  points., 
each  vanish  identically.  This  important  condition  must  hold  for  a solution 
to  exist  to  the  linear  algebraic  system  for  the  discretized  pressure  field; 
it  is  discussed  further  in  the  next  section  (see  Eqs.  (29)  and  (30)). 
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IV.  THE  SOLUTION  OF  THE  PRESSURE  EQUATION 

At  each  time  step  in  the  fully  nonlinear  problem  it  is  necessary 
to  calculate  the  solution  to  the  non-separable  elliptic  equation 

on  the  rectangle  R,  subject  to  the  Neumann  boundary  conditions 

= ?(*■>$) 

Here,  the  density  £>,  the  forcing  function  f , and  the  boundary  function  g 
all  depend  on  time,  but  are  known  to  us.  (in  the  linear  computation 
reported  here,  <p(x,  y)  is  replaced  by ^po(y),  which  yields  a separable 
elliptic  equation  as  noted  in  Section  II.  However,  the  solution 
procedure  is  developed  for  the  fully  nonlinear  problem. ) In  the 
discretization,  Equations  (25)  are  replaced  by  finite  difference 
equations,  Eqs.  ( l^d)  or  (2ld)  and  (23)  in  Section  3.  The  finite 
difference  equations  (Hd)  or  (5tJd)  then  directly  incorporate  the 
boundary  conditions  (23).  Hence,  the  calculation  of  the  pressure 
requires  only  the  solution  of  the  linear  equations  represented  by 
these  combined  equations. 

Physically  Equations  (25)  produce  a unique  solution  only  when  the 
pressure  is  prescribed  at  some  particular  point.  Equivalently,  there  is 
a unique  solution  with  mean  pressure  zero.  Given  any  particular  solution 
we  can  generate  a new  solution  by  adding  an  arbitrary  constant  to  the 
pressure.  The  linear  equations  (l9d),  (21d)  behave  analogously;  the 
linear  system  has  a family  of  solutions  such  that  there  is  a unique 
solution  with  mean  pressure  zero,  and  the  family  of  solutions  can  be 
generated  by  adding  constants  to  any  particular  solution. 

The  calculation  of  the  pressure  requires  us  to  successfully  meet 
the  constraints  of  the  larger  system  of  equations: 


(a. 5 o.) 
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the  solution  method  must  take  into  account  non- uniqueness ; the  solution 
method  must  he  able  to  solve  large  linear  systems  accurately,  since 
there  are  NM  equations  (l9d);  and  it  is  very  important  that  the  s lution 
be  obtained  quickly  since  the  calculation  is  made  at  each  time  step. 

The  finite  difference  equations  (l9d)  can  be  conveniently  represented 
in  matrix  notation  as 


Ap  = z 


where  p (p13_,  P^,  P-^,  ...  P^,  P^,  P22 , ...  P2M>  ...  ...  Pffl) 


is  a vector  of  dimension  NM,  as  it  z,  the  vector  of  right  sides  from  Eqs . 

(l9d)  or  (21d)  (modified  by  boundary  conditions  (23)).  The  matrix  A gives 

the  coefficients  from  the  density  and  the  finite  differences.  Here  the 

subscript  T denotes  transpose.  Although  A is  of  dimension  NM  by  NM,  it 

has  at  most  5 non-zero  elements  in  any  row. 

As  noted  above,  the  solution  to  Eqs.  (25)  is  determined  only  to  within 

an  additive  constant.  Similarly,  the  solution  to  the  system  of  linear 

equations  (2o)  is  also  determined  only  to  within  an  additive  constant. 

The  non-uniqueness  of  the  solution  is  represented  by 

A (p  + ae)  = z 
T 

where  e = (l,  1,  1,  1 ...  l)  ; a is  any  scalar  and  p any  solution  of  (26). 
This  can  be  represented  also  as : 

| Ae  = 0. 

| which  is  the  statement  that  e is  an  eigenvector  of  A corresponding  to  zero 
eigenvalue.  The  solution  p with  mean  pressure  zero  is  the  solution  of  (26) 
which  satisfies 


-T 

p e = 0. 


The  existence  of  a solution  requires  the  following  equation  to  hold 


z e = 0. 


(26) 


(27) 


(28) 


(29) 


(30) 


20 


This  condition  on  the  difference  equations  (l9d)  or  (21d) , incorporating 
boundary  conditions  (23b),  is  simply  the  statement  that  the  right  hand 
sides  sum  to  zero.  As  noted  at  the  end  of  the  last  section,  see  particularly 
Eq.  (24)  and  the  surrounding  discussion,  this  condition  is  guaranteed  by  the 
formulation  of  the  linearized  difference  equations.  It  also  must  be 
guaranteed  when  the  nonlinear  difference  equations  are  developed.  We  assume 
equation  (30)  and  proceed  to  calculate  the  solution  to  the  linear  system 
(26)  which  satisfies  the  mean  pressure  requirement  (29). 

Standard  direct  factorization  methods  (Gaussian  elimination)  for 

2 

calculating  the  solution  to  (2 6)  are  infeasible,  because  0 ((NM)  M) 

operations  at  best  are  necessary.  Iterative  methods  are  often  used  to 

solve  sparse  linear  systems:  the  method  of  successive  over- relaxation 

(SOR)  is  probably  best  known.  However,  a successful  implementation  of  SOR 

requires  at  least  a good  estimate  of  the  over- re taxation  parameter.  The 

convergence  rate  even  with  an  optimal  parameter  still  yields  0 ((NM)N  log  N) 

22 

operations  for  a solution  of  (26). 

The  specific  structure  of  the  finite  difference  operation  on  the 
rectangle  R makes  possible  special  fast  direct  methods  for  the  solution 
of  separable  elliptic  equations.  The  algorithms  discussed  in  Schwartztrauber 
and  Sweet"*''1'  are  the  most  current  extensions  of  the  ideas  behind  the  fast 
finite  Fourier  transform.  However,  even  these  general  algorithms  require 
separability  (the  coefficients  of  the  derivatives  with  respect  to  x involve 
only  x;  the  coefficients  of  the  derivatives  with  respect  to  y are  functions 
only  of  y),  and  the  solution  procedure  which  will  be  presented  here  is  for 
the  case  of  interest  in  nonlinear  problem.  The  pressure  equation  can  be 
transformed  to  separable  form  only  if  is  known  analytically  and  has  two 
continuous  derivatives.  Since  ^ is  known  only  on  the  staggered  grid,  these 
conditions  are  not  met. 
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The  solution  method  we  have  adopted  is  a hybrid  method  which  combines 
an  iteractive  algorithm,  conjugate  gradients  (denoted  c.g.),  with  a fast 
direct  Poisson  solver.  The  conjugate  gradients  algorithm  provides  an 
iterative  technique  for  solving  Ap  = z , provided  A is  a positive  definite 
symmetric  matrix  (equivalently,  all  eigenvalues  of  A are  positive).  Each 
iteration  of  conjugate  gradients  requires  one  matrix-vector  product,  A 
times  x or  roughly  5 MN  operations.  However,  the  convergence  rate  depends 
critically  on  the  ratio,  k = X^/X^,  the  largest  eigenvalue  divided  by 
the  smallest.  The  ratio  is  large  for  the  linear  system  (26),  and  con- 
vergence is  fast  only  if  k is  close  to  unity.  To  improve  (drastically) 

12 

1 the  rate  of  convergence,  we  write,  following  Concus  and  Golub, 

A = A + A 
s n 

where  A and  A are  matrices  chosen  so  that  A ^z  can  be  calculated  quickly 
s n s 

with  a fast  Poisson  solver  and  the  entries  in  A are  small  compared  to  A . 

n s 

The  best  choice  of  Ag  is  not  known  theoretically.  Our  solution  is 
to  take 

A = D LD 
s 

where  L is  the  finite  difference  matrix  using  a five-point  difference 

1 2 

formula  for  the  Poisson  equation  V P = f with  Newmann  boundary  conditions. 

Here  D is  a positive  diagonal  matrix  chosen  so  that  diagonal  (A  ) = diagonal 

s 

(A).  The  remainder  A^  is  a matrix  with  zero  main-di agonal , and  non-zero 
entries  present  on  only  four  off-diagonals.  Such  a choice  is  justified 
on  empirical  grounds  only  if  the  conjugate  gradient  convergence  rate  is 

I 

sufficiently  fast.  For  our  test  problems  this  has  certainly  been  the  case. 
There  is  also  some  theoretical  basis  for  this  choice;  see  Concus  & Golub ^ . 
Then  we  formally  replace  the  equation 

Ap  = z 
by 


(D  AD)  a = Dz 


(31) 


(34) 


(32) 
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where  DAD  is  a new  "scaled"  matrix  and  Dz  is  the  corresponding  scaled 
right  hand  side.  We  then  solve  the  scaled  problem  using  conjugate  gradients. 
The  solution  p of  (26)  is  then  Da. 

If  we  have  made  a suitable  choice  for  A , the  eigenvalue  ratio  X /X 
for  the  matrix  D AD  will  be  very  close  to  the  optimum  value,  unity.  The 
solution  of  (32)  requires  one  matrix- vector  product  Ax,  and  the  solution 

of 

As  x = y 

at  each  iteration.  (The  matrix  D does  not  appear  explicitly  in  the  actual 

algorithm. ) The  details  of  this  general  approach  are  given  in  Concus  , 

13 

Golub  and  O'Leary.  Although  the  solution  of  (33)  requires  0 (NM  log  N) 
operations,  the  entire  process  will  be  efficient  if  we  need  to  solve  (33) 
only  a very  few  times.  Typical  calculations  for  our  code  require  no  more 
than  three  to  five  Poisson  solutions  to  obtain  the  solution  of  (26). 

The  only  difficulty  with  the  approach  outlined  above  is  the  non- 
uniqueness of  the  solution.  The  matrix  A which  represents  the  finite 
difference  equations  is  only  positive-semi definite , since  e is  an 
eigenvector  corresponding  to  the  eigenvalue  0.  However,  all  other 
eigenvalues  are  positive.  The  iterations  in  the  conjugate  gradients 
algorithm  will  generate  a solution  orthogonal  to  e if  and  only  if  the 
iterations  are  begun  with  a vector  orthogonal  to  e. 

The  natural  solution  produced  by  the  method  is  the  zero  mean  pressure 
solution;  any  other  solution  can  be  produced  from  this  particular  solution 
by  adding  a constant. 

There  is  an  additional  complication  in  the  hybrid  algorithm  since 
both  A and  L,  the  discrete  Poisson  operator,  are  singular;  Ae  = Le  = 0. 

The  iteration  matrix  D AD  has  a single  zero  eigenvalue  with  eigenvector 
D-1e.  In  we  show  that  for  a more  general  splitting  of  matrix  A, 
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the  iteration  matrix  corresponding  to  D AD  has  the  desired  single  zero 

eigenvalue  with  both  A and  the  matrix  corresponding  to  D singular  if  and 

only  if  the  corresponding  eigenvectors  of  A and  the  more  general  D are  not 

T -1  1 

orthogonal.  This  can  he  assumed  in  our  calculation  since  e D e = — — 
i 

+ - — + . . . + — — and  all  the  d.  . are  positive. 
d12  dNM 

It  is  shown  in  [l4]  that  a choice  of  Ag  with  the  same  eigenvector  of  zero 
eigenvalue  as  A is  optimal  in  some  respects,  hut  such  a choice  seems 
difficult  to  realize.  Our  simpler  choice,  as  outlined  above,  proves 
entirely  satisfactory.  However,  we  must  note  that  the  search  direction 
in  the  c.g.  iteration  must  now  he  orthogonal  to  D ^e,  not  e.  This  requires 
that  the  solution  to  the  discrete  Poisson  equation  Lx  = y he  normalized 
to  mean  pressure  zero.  This  normalization  must  he  done  by  the  user  of  the 
Schwarztrauher-Sweet  routines^1  and  probably  is  also  necessary  for  other 


fast  Poisson  solvers. 


24 


V.  Computational  Procedure  for  Solving  Linear  Finite  Difference  Equations 

In  Section  III  the  finite  difference  equations  describing  linearized 
heating,  Eqs.  (I9)-(23),  were  presented,  and  in  Section  IV  the  procedure 
used  to  solve  the  pressure  equation  was  described.  In  this  section  the  solu- 
tion procedure  for  the  complete  set  of  linearized  difference  equations  will 
be  presented  and  some  aspects  of  the  computation  will  be  discussed. 

Since  the  finite  difference  scheme,  Eqs.  (19)  is  a three  level  scheme, 
values  of  velocities  and  density  at  two  levels  of  time  are  required  to 
initiate  the  computation.  Equations  (21),  a scheme  first  order  in  time. 


p _ given  initial  data.  Equations  (.21;  are  solved  in  the  following  order. 
First  the  densities  p^.  are  obtained  from  Eq.  (2la) . Then  the  pressures 


u_^_. , v_  are  computed  from  (21b)  and  (2lc). 

One  point  in  this  process  is  noteworthy.  The  initial  pressure  data. 


pressure  p are  obtained  using  p_.  Such  a procedure  was  adopted  because 
our  analytical  studies  of  the  finite  difference  scheme  first  order  in  time 
demonstrated  that  it  was  most  appropriate  for  reproducing  the  test  problem 
(see  Reference  10  for  details) . 

In  subsequent  time  steps  the  order  in  which  the  difference  equations 


ij 

p?.  are  obtained,  using  the  p^.,  from  Eq.  (21d).  Finally  the  velocities 


1 1 


provided  to  the  computational  scheme  is  used  only  as  an  initial  "guess"  for 


the  iterative  pressure-solving 


0 


procedure:  the  initial  values  of  the 

1 „ , , , _ , t 


are  solved  is  changed  somewhat.  First,  the  pressure  equations,  Eq.  (l9e) 
for  the  mean  pressure  and  Eq.  (19 d)  for  the  pressure  variation,  are  solved. 
Then  the  density  and  the  two  velocities  are  updated  using  Eqs.  (19a)- (19c). 
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It  should  be  noted  that  in  both  the  first  time  step  and  subsequent 
ones,  the  solution  to  the  pressure  equation  assures  that  the  velocity  com- 
ponents will  satisfy  the  proper  divergence  condition.  Equation  (5)  is  the 
imposed  condition  on  the  divergence  of  the  velocity  in  the  continuous  case; 
this  equation  also  applies  in  the  linear  problem,  with  the  dimensionless 
D(x,y,t)  being  specified  by  the  fourth  of  Eqs.  (9).  The  discrete  version 
of  Eq.  (5)  is 


1 f n+1  n+1  x 1 , n+1  n+1  x _n+l 

— (u.  . - u.  .)+  — (v..  -v.  . .)  = D 

fix  i,j  x-1, j 6y  ij  1,3-1  il 


(35) 


where  ]>”  is  given  by  Eqs.  (2 0 ) * This  condition  is  very  important  since 
it  specifies,  in  the  discrete  (and  linear)  case,  what  effect  the  heat  source 
imposes  through  continuity  upon  the  velocity  field.  Substituting  for  u?*\ 
v?T^  from  Eqs.  (19b)  and  (19c),  and  using  Eq.  (35)  again, with  n-1  replacing 
n+1,  yields  the  equation  for  the  pressure,  Eq.  (25d)  for  all  steps  after  the 
first.  Hence  proper  solution  of  this  pressure  equation  assures  that  the 
velocities  at  time  step  n+1  will  satify  the  proper  divergence  condition, 

Eq.  (35).  An  analogous  argument  for  the  first  time  step,  using  Eq.  (2lb), 

(21  ) and  (35),  shows  that  the  pressure  equation  (21d)  assures  the  proper 
divergence  condition  at  the  first  time  step. 

When  there  is  no  imposed  heating,  so  that  = 0,  Eqs.  (19) 

and  (21)  are  homogeneous  (as  noted  earlier).  In  Reference  10  the  homogeneous 
equations  (19)  with  the  first-order  starting  procedure  implied  by  the  homo- 
geneous form  of  Eqs.  (21)  were  solved  analytically  for  internal  waves  in  a 
rectangular  enclosure.  Comparisons  were  made  between  this  analytical  solu- 
tion to  the  finite  difference  equations  and  one  generated  by  computational 
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procedures.  These  comparisons  showed  that  the  accuracy  with  which  the 
difference  equations  were  solved  computationally  was  determined  by  the 
accuracy  specified  for  the  solution  of  the  pressure  equations  (19d)  (or 
(2ld)  on  the  first  time  step) . When  the  accuracy  specified  for  the  pressure 
equation  was  10  the  analytically  and  computationally  generated  solutions 
agreed  to  1 or  2 parts  in  10^. 

Further  comparisons  were  made  between  the  analytical  solution  to  the 
finite  difference  equations  and  the  analytical  solution  to  the  continuous 
problem,  the  homogeneous  form  of  Eqs.  (13).  These  solutions  differ  by  the 
truncation  error  resulting  from  replacing  differential  equations  by  differ- 
ence equations.  The  truncation  errors  in  turn  depend  upon  the  spatial  grid 
size  relative  to  characteristic  spatial  variations  in  the  solution  and  upon 
the  time  step  size  relative  to  characteristic  temporal  variations  in  the 
solution.  Internal  waves  in  a rectangular  enclosure  were  examined  analytically 
to  test  the  finite  difference  schemes  and  computational  procedure.  Such 
waves  are  expected  once  the  fire  has  stratified  the  gas  in  the  room. 

A single  internal-wave  mode  was  specified  initially  (in  either 
continuous  or  discretized  form)  and  its  time  history  was  calculated  by  way 
of  the  differential  or  difference  equations.  The  mode  was  specified  by 
two  integers,  the  number  of  half  wave  lengths  in  the  x-direction  and  the 
number  of  half  wavelengths  in  the  y-direction  filling  the  rectangular 
enclosure  and  satisfying  boundary  conditions.  Specification  of  the  mode 
then  determined  the  eigenf requency , through  the  equations  of  motion,  with 
which  this  characteristic  mode  oscillated.  The  stability  of  the  finite 
difference  equations  was  found  to  require  that 
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This  condition  requires  that  the  time  step  size  be  a fraction  of  the  Brunt- 
Vaisala  period  for  the  ambient  fluid,  eigenf requencies  being  inversely 
proportional  to  this  period.  Accuracy  of  the  finite  difference  equations 
was  found  to  require  that  the  mode  number  relative  to  the  mesh-cell 
number  M in  the  x-direction  be  small,  the  mode  number  relative  to  the 
mesh-cell  number  N in  the  y-direction  be  small  and  the  time  step  size 
relative  to  the  Brunt-Vaisala  period  be  small. 

As  an  example,  a 5 by  2 mode  internal  wave  was  calculated  for  an 
ambient  stratification  parameter  (see  equation  16a)  equal  to  unity. 

For  this  stratification,  the  Brunt-Vaisala  period  is  2tt  in  dimensionless 
time  units.  A 31  x 30  mesh  was  employed,  with  6t  = 0.25.  At  a dimensionless 
time  t=10,  the  discrete  and  continuous  solutions  differed  by  about  3%. 

During  this  time  the  discrete  solution  oscillated  about  the  continuous 
result  in  an  irregular  fashion.  This  behavior  in  the  discretization  error 
is  characteristic  of  the  spurious  computational  mode  introduced  by  the 
leapfrog  time  differencing  procedure.  The  continuous  eigenvalue  a which 
determines  the  actual  period  of  the  mode  under  consideration  is  a=. 92807 
for  this  example.  The  corresponding  discrete  mode  eigenvalue  o*  for  this 
problem  may  be  analytically  calculated'*’ to  be  a*-. 923955.  Thus,  the 
computational  scheme  should  be  able  to  follow  oscillations  sufficiently 
accurately  to  avoid  errors  due  to  phase  shifts  between  modes  for  times 
t 0(100).  Throughout  the  time  interval  0 £ t < 10,  the  finite  difference 
equations  were  solved  computationally  to  a few  parts  in  10^.  This  estimate 
is  based  on  direct  comparison  between  the  analytically  obtained  and  numerically 
computed  results. 
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Last,  we  mention  briefly  computational  times  required.  All  computa- 
tions reported  in  the  next  section  (and  those  performed  for  results 
reported  in  Reference  2)  were  carried  out  on  the  N.B.S.  UNIVAC  1108.  Most 
computations  were  performed  using  M = 31  and  N = 31  for  the  mesh  and  were 
run  for  up  to  two  hundred  time  steps.  For  a calculation  performed 
starting  with  a stratified  fluid,  between  three  and  four  iterations  on 
average  were  taken  by  the  conjugate  gradients  portion  of  the  pressure 
solver  and  the  computer  time  required  per  time  step  was  between  4 and  5 
seconds.  For  example,  one  computation  of  70  time  steps  took  302  seconds 
of  CPU  time  and  another  took  3l6  seconds.  For  a calculation  in  an 
essentially  unstratified  fluid  (in  which  the  elliptic  equation  for  the 
pressure  becomes  Poisson's  equation),  only  one  iteration  was  required  and 
about  2 seconds  per  time  step  was  taken.  In  this  case,  a computation  run 
for  200  time  steps  required  402  seconds  of  CPU  time. 

Computations  on  smaller  meshes  required  substantially  less  time.  For 
example,  a computation  performed  for  a stratified  fluid  using  M = 7 and  N = 8 
required  66  seconds  of  CPU  time  to  run  200  time  steps.  Another  calculation, 
performed  for  an  unstratified  fluid  with  M = 15  and  N = l6  and  also  run 
for  200  time  steps  required  122  seconds. 

The  computational  times  quoted  do  not  include  CPU  time  required  to 
produce  the  computer-generated  plots  discussed  in  the  next  section.  When 
plots  were  generated,  the  results  were  read  to  file  during  the  computation 
and  processed  later  to  produce  the  graphics.  Computer  time  required  for  , 
the  graphics  will  be  discussed  in  the  next  section. 
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VI.  Computer-Generated  Results 

Because  Egs.  (13)  with  boundary  conditions  (14)  and  (15)  are 
linear,  they  can  be  solved  analytically  (and  these  expressions  evaluated 
numerically)  or  by  a combination  of  analytical  and  computational  tech- 
niques. Unfortunately,  neither  of  these  procedures  has  been  implemented 
as  a computer  code  at  this  time  so  that  comparison  of  the  results  reported 
in  this  section  with  results  obtained  by  an  alternate  procedure  cannot  be 
made  yet.  Presently,  we  are  implementing  one  approach  using  a mixture  of 
analytical  and  numerical  techniques,  and  we  will  make  such  comparisons 
later . 

A.  Graphics 

From  the  finite  difference  equations  (19),  we  see  density,  horizontal 
and  vertical  velocity  components  and  pressure  are  all  computed  dependent 
variables  at  each  time  step.  Typically  in  the  computations  we  have  run, 
we  have  used  a grid  composed  of  31  x 31  mesh  boxes:  hence  at  each  time 

we  have  4 dependent  variables,  each  being  calculated  at  approximetely  1000 
points.  Clearly  with  such  a large  number  of  values  computed,  graphical 
display  is  required  to  digest  the  results.  Generally  we  have  used  two 
computer  plotting  routines  to  allow  us  to  visualize  results;  both  routines 
were  developed  by  the  Mathematical  Analysis  Division  specifically  for  this 
application  and  have  been  found  to  be  essential  for  rapid  assimilation  of 
results.  Most  plots  shown  in  this  report  are  direct  reproductions  of  these 
computer-generated  plots. 

Both  plotting  routines  display  the  dependent  variables  on  the  rectangular 
region  of  computation.  One  routine  displays  profiles  of  the  horizontal  and 
vertical  velocity  components.  The  other  plotting  routine  displays  contours 
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of  constant  density  (isopycnic  contours). 

The  horizontal  velocity  is  plotted  as  a function  of  vertical  position 
at  prescribed  horizontal  locations,  together  with  plots  of  vertical 
velocity  as  a function  of  horizontal  position  at  prescribed  vertical 
locations.  (The  locations  at  which  horizontal  or  vertical  velocity  plots 
are  desired  can  be  specified.)  These  superimposed  plots  give  a clear 
picture  of  the  velocity  field  at  any  instant  of  time. 

At  a horizontal  location,  this  plotting  routine  takes  the  values  of 
the  computed  horizontal  velocity  components  at  the  vertical  mesh  points 
and  fits  this  data  with  cubic  splines  using  a package  "splines  under 
tension"  written  by  A.K.  ClineV  * iince  the  vertical  mesh  points  at  which 
the  horizontal  velocity  components  are  computed  is  staggered  with  respect 
to  the  grid  lines  covering  the  computational  region,  horizontal  velocities 
are  not  computed  at  the  top  or  at  the  bottom  boundaries.  The  last  computed 
horizontal  velocity  is  one  half  mesh  point  from  the  boundary,  and  the 
velocity  profiles  are  extrapolated  to  the  boundaries  using  the  fitted  cubic 
splines  (not  under  tension) . 

Similarly  at  the  vertical  location,  the  plotting  routine  fits  the 
computed  vertical  velocity  values  as  a function  of  horizontal  position  with 
cubic  splines.  However,  for  the  vertical  velocity,  the  normal  derivative 
at  the  left  and  right  walls  can  be  shown  to  be  zero  for  the  linearized 
equations . 

Using  this  condition  and  the  values  of  the  vertical  velocity  at  the 
two  points  adjacent  to  each  wall,  a quadratic  equation  was  used  to  extrap- 
olate the  value  of  the  vertical  velocity  at  each  wall.  Then  the  package 
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written  by  A.K.  Cline  was  used  to  fit  the  vertical  velocity  data  and  the 
end  point  information  by  cubic  splines  (not  under  tension) . An  earlier 
version  of  this  package  is  described  in  references  15  and  16. 

Finally,  the  velocity  plots  shown  below  were  prepared  with  an  inter- 
active Fortran  program  using  Tektronix  Plot  10  Advanced  Graphing  II  soft- 
ware and  Terminal  Control  System  (TCS)  software^’ Hie  display  device  is  a 
Tektronix  4014  display  terminal. 

The  second  plotting  routine  gives  constant-density  contours.  It  first 

creates  an  intermediate  code  with  an  interactive  Fortran  using  the  NCAR 
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graphics  software  CONREC  and  DASHCHAR  together  with  the  system  plot  package 

and  the  Tektronix  Plot  10  TCS.  Then  the  intermediate  form  is  translated 
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using  the  NCAR  portable  graphics  translator  written  by  T.  Wright  and  the 
Tektronix  Plot  10  TCS  software.  The  plots  are  displayed  on  a Tektronix 
40l4  display  terminal.  Typical  times  required  to  generate  a velocity  or 
a density  plot  are  five  to  fifteen  seconds  of  CPU  time. 

B.  Linear  Heating  in  an  Ambient  Homogeneous  Fluid 

When  the  ambient  fluid  is  homogeneous,  (y)  is  a constant  and  Eqs. 
(13)  simplify  considerably.  For  this  case  it  is  valuable  to  examine  the 
analytical  form  of  the  solution  even  though  the  solution  is  not  calculated 
for  comparison  with  the  numerical  results.  Several  observations  about  the 
nature  of  the  flow  field  can  be  made  from  the  analytical  solution,  and 
these  observations  are  useful  discussing  the  numerical  results. 

When  pQ  is  constant,  the  equation  for  the  density  in  Eqs.  (13)  can  be 
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solved  explicitly  using  the  separable  form  for  Q (x,y,t)  given  in  Eq.  (17). 

Acj.]  - flU  A 


where 


t>*Uvd  ''  Sr'A  vrr 


and  where  the  density  perturbation  from  is  initially  taken  to  be  zero. 
The  density  is  found  to  be  a function  of  time  multiplying  a function  of 
space,  and  density  can  be  expected  to  display  exactly  the  same  spatial 
behavior  at  different  times  with  only  the  amplitude  changing  from  one 
time  to  the  next. 

The  velocity  field  can  be  decomposed  into  two  fields,  one  derivable 
from  a potential  the  other  from  a stream  function: 


u. 


= £>Al 

"&x 


’ 


P'V' 


where  <p  is  the  potential  and  ip  stream  function.  The  boundary  conditions 
on  these  two  functions  arise  from  the  condition  that  the  normal  velocity 
at  a boundary  be  zero  (no  outflow  at  a boundary),  Eqs . (10):  they  are 

^ i ■ t ^ 

— = 0 and  i|)  = 0 on  the  boundary,  x = 0,  and  y = 0,  1 (where  — means 

the  derivative  normal  to  the  boundary). 

The  continuity  equation,  Eq . (5),  implies 

va+  = h (x.,#.,  t)  = fft)  b*  (*,#) 

If  we  denote  by  <£  (*  ^ y J the  solution  to  the  equation 

V7 - D*(x>y) 

with  — O at  and  ^ ~Oi  i , then 
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We  define  the  vorticity  to  be 

C v =r 


a>r 


Then 


The  equation  for  the  time  evolution  of  the  vorticity  is  obtained  by 
taking  the  curl  of  the  momentum  equations  in  Eqs . (13) 

1W  - -J-  \L 

The  solution  is 

3 X ^ 7 

The  initial  conditions  that  u (x,  y,  0)  = v (x,  y,  0)  = 0 have  been  used 
to  determine  that  the  initial  vorticity  distribution  is  zero. 

If,  as  before,  we  denote  by  J (x,  y)  the  solution  to  the  equation 

with  | = 0 at  x = 0,  ^^5  and  y = 0,  1,  then 

The  point  to  note  from  the  analytical  forms  derived  above  is  that 
the  temporal  behavior  of  the  velocity  field  obtained  from  the  potential 
is  different  from  that  obtained  from  the  stream  function.  In  particular, 
when  f (t)  grows  linearly, 


(-t)  oC  t 
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and 


Hence,  initially  the  potential  contribution  to  the  velocity  field  will 
dominate  and  the  velocity  field  will  grow  approximately  linearly  with 
time, whereas  later  the  contribution  from  the  stream  function,  or  rather 
due'  to  the  generation  of  vorticity,  will  dominate  and  the  velocities  will 
grow  roughly  as  the  cube  of  time. 

In  Figure  2-6  plots  of  horizontal  and  vertical  velocity  are  dis- 
played for  a case  of  heating  in  an  unstratified  fluid.  The  heat  source 
is  of  the  form  given  in  Eq  (17)  with 


Values  of  the  parameters  for  the  heat  source  shown  are  3 = 24.5,  X = 4, 
x^  = 0.5  (cf  Eq.  (17))  and  A=  0.1:  hence  the  source  is  centered  along 

the  bottom  of  the  floor,  its  strength  diminishing  more  rapidly  laterally 
than  vertically. 

Figure  2 shows,  at  time  t = 0.1,  horizontal  velocity  as  a function 
of  vertical  position  at  three  locations  (x  = 0.25,  0.48  and  0.74)  and 
vertical  velocity  as  function  of  horizontal  position  at  three  locations 
(y  = 0.25,  0.48,  0.74).  (Note  that  the  heat  source  and  the  flow  field 
are  symmetric  about  the  centerline,  but  the  horizontal  positions  used 
for  plotting  are  not. ) Figure  3 shows  the  same  horizontal  and  vertical 
velocity  plots  at  t = 0.2.  Comparison  of  these  two  figures  shows  that 
the  velocities  in  Fig.  3 have  the  same  spatial  form,  but  are  nearly 
double  those  in  Fig.  2.  At  "early  time"  during  the  linear  calculation, 
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when  these  two  plots  were  made,  the  velocities  are  dominated  by  the 
potential  flow  as  the  analysis  above  indicates.  The  potential  flow 
represents  an  expansion  of  the  gas  due  to  heating  by  the  source.  At  this 
time  the  velocities  grow  linearly  with  time,  and  Figs.  2 and  3 confirm 


this  behavior. 
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Figures  4 and  5 show  these  horizontal  and  vertical  velocity  plots 
at  t = 2.0  and  t = 4.0  respectively.  Note  the  change  in  the  velocity 
scales  from  Figs.  2 and  3.  At  these  times,  which  are  "late  time"  in 
the  linear  flow,  the  flow  field  is  dominated  by  the  vorticity  field  for 
which  the  physical  significance  is  discussed  below.  The  velocity  plots 
in  Fig.  5 have  the  same  spatial  form  but  are  approximately  eight  times 
larger  than  those  given  in  Fig.  4,  a fact  consistent  with  the  analysis 
above,  which  shows  that  the  stream  function  (produced  by  the  vorticity 
field)  grows  as  the  cube  of  time. 

are  contributions  to  the  velocities  from  both  the  potential  and  the  stream 
function;  hence  t = 0.7  is  an  "intermediate  time"  for  the  linear  calculation 
in  which  the  velocity  field  is  neither  dominated  by  the  potential  nor  by 
the  vorticity.  (Note  that  Fig.  6 has  the  same  scales  for  velocity  as  Figs. 

2 and  3 . ) 

An  observation  should  be  made  at  this  point  regarding  the  qualitative 
behavior  of  the  linear  flow  field  exhibited  in  Figs.  2-6.  The  overall 
nature  of  the  flow  field  is  what  one  would  expect  (and,  for  the  later 
times,  qualitatively  what  people  find  in  a room  fire).  During  the  "early 
time"  of  the  linear  calculation  as  noted  before  (when  the  velocities  are 
dominated  by  the  potential  function),  the  flow  is  strictly  outward  from 
the  heated  region  as  the  gas  is  heated  and  expands.  At  "later  time" 
during  the  linear  calculation,  the  flow  is  upwardly  directed  above  the 
heated  region  and,  by  continuity,  down  along  the  sides  of  the  enclosure. 
Qualitatively,  this  is  the  behavior  observed  in  an  enclosure  fire.  In 
addition,  in  the  region  above  the  most  intensive  heating  (approximately 
y £ 0.17  here),  the  horizontal  flow  is  mildly  outwardly  directed  from 
the  region  above  the  heat  source  (the  plume  region).  Near  the  bottom  of 
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the  enclosure,  however,  there  is  a relatively  strong  inflow.  In  room 
fires  a strong  inflow  or  "floor  jet"  is  often  observed  near  the  floor. 

Such  a qualitatively  good  picture  from  a linear  calculation  is  very 
encouraging. 

Finally,  in  Figs.  7 and  8 plots  of  constant- density  contours  are 
shown  at  two  different  times,  t = 0.5  and  t = 2.0  respectively.  Two  points 
should  be  noted.  First,  as  observed  above  from  the  analytical  solution, 
the  spatial  behavior  of  the  density  at  the  two  different  times  is  the  same, 
only  the  magnitude  of  the  density  perturbation  has  changed  between  the  two 
times.  Second,  we  note  that  the  spatial  behavior  of  the  density  perturbation 
is  proportional  to  D (x,  y)  and  hence  directly  reflects  the  nature  of  the 
specified  heat  source  Q (x,  y) . 
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C.  Linear  Heating  in  a Stratified  Fluid 

When  the  ambient  fluid  is  stratified  with  an  exponential  variation 
as  given  in  Eq.  (16a),  a formal  solution  to  Eqs . (13)  can  still  be  obtained 
(because  the  equations  are  linear) . However,  because  of  the  complexity  of 
this  solution,  it  is  not  very  helpful  for  gaining  an  understanding  of  the 
behavior  of  the  flow  field.  Therefore,  in  this  section  results  of  some 
computations  solving  the  complete  difference  equations,  Eqs.  (19)  and  (21), 
will  be  presented  and  discussed,  but  no  additional  analysis  will  be  presented. 

In  all  of  the  computations  reported  in  this  section,  a reasonably 
refined  mesh,  I = J = 31,  was  chosen.  The  enclosure  was  taken  to  be  a 
square,  /'i  = 1.0,  and  the  gas  constant  = 1.4.  Parameters  describing  the 
heat  source  were  6 = 24.5,  \ = 4.0  in  Eqs.  (19)  and  A = 0.1  in  Eqs.  (20). 

The  parameter  e,  specifying  the  relative  error  to  which  the  pressure  equation 
is  to  be  solved,  was  taken  to  be  10  ^ , and  the  dimensionless  time  step  was 
taken  to  be  6t  = 0.2.  All  other  parameters  relevant  to  the  individual 
computations  reported  will  be  stated  when  the  results  are  discussed. 

An  interesting  observation  can  be  made  at  this  point  regarding  the 
relative  effort  required  to  obtain  numerical  solutions  (in  the  form  of  plots 
say)  as  the  complexity  of  the  (linear)  problem  increases.  In  a previous  paper, 
Ref.  10,  the  authors  obtained  solutions  to  the  equations  for  internal  waves 
in  enclosures.  In  the  continuous  case,  these  waves  are  solutions  to  the 
homogeneous  form  of  Eqs.  (13),  and  in  the  discrete  case,  they  are  solutions 
to  the  homogeneous  forms  of  Eqs.  (19)  and  (21).  Analytical  solutions  to 
both  the  continuous  and  the  discrete  cases  were  obtained,  and  a numerical 
solution  was  also  obtained  by  the  time-stepping  procedure  described  in  this 
report.  To  generate  a code  and  compute  numerical  values  for  the  homogeneous 
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problem  was  very  simple  using  the  analytical  solutions,  but  required  much 
effort  for  the  direct  numerical  time-stepping  procedure.  Once  the  time- 
stepping  procedure  was  implemented  as  a code,  it  was  a, easy  task  to  modify 
it  to  calculate  the  inhomogeneous  problems  with  constant  or  variable  ambient 
density  (described  in  the  last  section  and  this  section  respectively).  To 
calculate  values  from  the  analytical  solution  to  the  inhomogeneous  problem 
with  constant  background  density  in  the  continuous  case  requires  some 
effort,  but  not  too  much.  However,  to  perform  a similar  task  when  the 
background  density  is  variable  is  a rather  big  effort.  In  short,  as  the 
complexity  of  the  problem,  even  though  linear,  increases,  the  effort  re- 
quired to  obtain  and  evaluate  analytical  results  may  become  much  greater 
than  straightforward  numerical  integration  by  difference  methods.  In  the 
nonlinear  case,  of  course,  analytical  solution  will  not  even  be  possible. 

Cil.  Heating  in  a Stratified  Fluid  by  a Centered  Heat  Source 

Internal  waves  can  be  sustained  in  a stratified  fluid.  As  noted 
above  these  waves  are  solutions  to  the  homogeneous  form  of  Eqs.  (13).  The 
heating  source  is  the  inhomogeneous,  driving  term  for  these  equations,  and 
this  source  excites  internal  waves  during  heating.  One  of  the  most  inter- 
esting features  of  the  plots  shown  in  this  and  the  subsequent  two  examples 
in  this  section  are  the  internal  waves  generated  by  the  heat  source.  With 
each  of  the  examples  other  observations  will  be  madp. 

In  this  first  example  the  stratification  parameter  Y = 1.0.  Hence 
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the  buoyancy  period  T = 2tt  {-g  ~ = 2tt,  and  to  observe  anY  effects 

of  internal  waves,  the  calculations  must  be  performed  for  at  least  a period 
or  two.  This  computation  was  run  for  14  dimensionless  time  units  with 
6t  = 0.2  so  that  70  time  steps  were  taken. 


In  Fig.  9 plots  of  horizontal  and  vertical  velocities  at  dimensionless 
time  t = 4.0  for  this  case  are  presented.  These  plots  should  be  compared 
with  those  presented  in  Fig.  5:  qualitatively,  the  plots  are  the  same. 

However,  some  quantitative  differences  can  be  observed.  First,  all  velocities 
are  somewhat  smaller  in  Fig.  9 (the  stratified  fluid)  than  in  Fig.  5;  the 
difference  is  particularly  noticable  in  the  vertical  velocities.  This 
effect  is  directly  traceable  to  the  stable  background  stratification  which 
tends  to  suppress  buoyant  vertical  velocities  induced  by  the  heating. 

Another  difference  appears  in  the  horizontal  velocities  where  the  transition 
between  outflow  from  the  center  of  the  heated  region  and  inflow  occurs  lower 
in  the  enclosure  due  to  the  background  stratification. 

The  next  three  figures,  Figs.  10,  11  and  12,  show  the  subsequent  devel- 
opment with  time  (t  =8.0,  12.0  and  14.0  respectively)  of  the  velocity  pro- 
files. Note  the  change  in  the  velocity  scales  from  Fig.  5.  It  appears  from 
these  plots  that  the  vertical  velocity  changes  little  in  magnitude  over  the 
time  period  represented.  However,  the  horizontal  velocities  change  sub- 
stantially during  this  period.  The  maximum  horizontal  velocity  increases 
with  time. and  the  transition  between  outflow  from  the  center  of  the  heat 
source  to  inflow  occurs  lower  in  the  enclosure  at  each  time.  Other 
distortions  of  the  profiles  as  functions  of  time  are  also  apparent,  and 
these  changes  are  a result  of  the  background  stratification. 

The  nature  of  the  internal-wave  effects  can  be  seen  more  clearly  in 
plots  of  the  density  perturbation  (from  the  background  exponentially  varying 
density)  at  times  t = 2.0,  4.0,  8.0  and  12.0  in  Figs.  13  - 16.  In  Fig.  13 
constant  perturbation-density  contours  at  t = 2.0  are  shown:  clearly  these 

are  qualitatively  very  similar  to  the  perturbation  density  contours  for  the 
unstratified  case  shown  in  Figs.  7 and  8.  However,  at  later  times  the 
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density  contours  become  qualitatively  quite  different,  the  most  dominant 
feature  at  any  single  time  being  that  these  contours  are  dramatically 
flattened.  Internal- wave  effects  become  apparent  by  comparing  plots  at 
different  times.  Such  comparison  shows  that  the  point  (or  points)  of 
highest  density,  denoted  by  the  letter  H,  sloshes  back  and  forth  within 
the  enclosure.  Furthermore,  examination  of  the  zero-perturbation  contour 
shows  undulations  with  increasing  time. 
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C.2.  Heating  in  a Stratified  Fluid  by  a Noncentered  Heat  Source 
In  this  calculation  the  center  of  the  heat  source  was  shifted  to 
x = 0.306  (compared  with  x = 0.5  in  previous  calculations).  All  other 
parameters  remained  the  same  as  in  previous  calculations,  and  the  com- 
putation was  run  for  70  time  steps,  or  to  t = 14.0. 

Figures  17,  18  and  19  show  velocity  profiles  at  t = 4.0,  8.0  and  12.0 
respectively.  (Note  the  change  of  scale ,of  the  plots  between  Fig.  17  and 
Fig.  18.)  The  most  distinguishing  feature  of  these  plots  (and  of  this 
calculation)  is  the  asymmetry  introduced  by  this  noncentered  heat  source. 

As  in  the  previous  case,  the  magnitude  of  the  vertical  velocities  change 
little,  but  the  shape  of  these  plots  change  somewhat  as  functions  of  time. 

The  horizontal  velocity  profiles  change  both  in  magnitude  and  in  shape  as 
functions  of  time.  The  locations  of  the  transitions  between  outflow  from 
the  center  of  the  source  to  inflow  change  as  functions  of  time,  occuring 
lower  in  the  enclosure  as  time  increases. 

Figures  20  - 23  show  constant-perturbation  density  contours  at  t = 2.0, 

4.0,  8.0  and  14.0  respectively.  At  t = 2.0,  these  contours  are  qualitative 
similar,  although  displaced  in  the  enclosure  because  of  the  noncentered  heat 
source,  to  those  shown  in  Fig.  13.  At  later  times  the  contours  become  flattened, 
and  undulations,  as  described  in  the  previous  case,  occur. 
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C.3.  Heating  in  a Stratified,  Two-Layer  Fluid  by  a Noncentered  Heat 
Source 

Q 

In  the  final  example  all  parameters  are  as  described  in  the 
previous  case  except  that  the  background  density  stratification  changes 
abruptly  at  y = 0.4,  the  fluid  being  of  uniform  density  below  this 
location.  Figure  24  shows  plots  of  constant-perturbation-density  at  t = 6.0 
in  this  calculation.  This  figure  is  shown  only  to  demonstrate  the  quali- 
tativelv  different  behavior  which  occurs  when  the  background  density  is 
composed  of  two  layers,  the  upper  one  of  variable  density  and  the  lower 
one  of  constant  density  (homogeneous).  In  the  lower  layer  the  profiles 
are  qualitatively  similar  to  what  would  be  expected  for  heating  in  a 
homogeneous  fluid.  The  structure  in  the  upper  layer  is  qualitatively 
quite  different,  showing  what  appears  to  be  much  detailed  structure 
superimposed  upon  some  mild  undulations.  The  detailed  structure  should 
not  be  seriously  considered  because  the  magnitude  of  these  variations  is 
small.  On  the  basis  of  previous  studies.  Reference  10,  the  discretization 
errors  expected  for  the  mesh  M = N = 31  is  of  the  order  of  a few  percent. 

(it  could  be  somewhat  larger  because  the  background  density  has  a discon- 
tinuity in  its  gradient  between  the  upper  and  lower  layers. ) The  magni- 
tude of  the  variations  of  the  detailed  structure  is  smaller  than  the 


expected  discretization  errors. 
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VII.  Summary 

4 

In  an  earlier  paper  the  equations  of  motion  for  thermally  driven 
buoyant  flows  were  derived.  These  approximate  equations  are  applicable 
to  nondissipative,  highly  nonadiabatic , buoyant  flows  of  a perfect  gas  in 
which  the  heat  is  added  slowly  by  a prescribed  source.  The  equations 
were  derived  by  assuming  that  the  time  scale  associated  with  increase 
of  the  heat  source  is  long  compared  with  the  transit  time  of  an  acoustic 
signal  across  the  spatial  extent  of  the  source.  Such  an  assumption  is 
appropriate  for  buoyant  flows  produced  by  room  fires.  The  resulting 
equations  admit  internal-wave  motions  while  filtering  out  high-frequency 
acoustic  waves.  They  are  a generalization  of  the  Boussinesq  equations 
valid  for  nearly  isobaric  flows  with  large  density  and  temperature  differences. 
This  combination  of  phenomena  has  not  been  extensively  examined  in  studies 
of  fluid  flows. 

In  the  present  paper  these  nonlinear  equations  were  recast  in  a form 

in  which  they  could  be  integrated  by  finite  difference  techniques.  In 
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three  dimensions  the  recast  equations  consist  of  three  evolution  equations, 
for  density  and  the  three  components  of  velocity,  along  with  a linear, 
nonseparable  elliptic  equation  for  the  spatially  variable  portion  of  the 
pressure.  (The  mean  pressure  as  a function  of  time  is  obtained  by  solving 
an  ordinary  differential  equation) . The  formulation  in  which  these  fluid 
quantities  are  used  as  dependent  variables  is  often  refered  to  as  a primitive 
variable  or  primitive  equation  formulation. 

The  equations  were  then  made  dimensionless,  specialized  to  two 
dimensions  and  linearized.  Finite  difference  approximations,  second 
order  accurate  in  both  space  and  time,  were  made  to  the  linearized  equations. 
The  dependent  variables  are  density,  the  horizontal  and  vertical  velocity 
components  and  pressure,  and  these  variables  are  defined  at  various  positions 
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within  a grid  cell  (see  Figure  lb) . 

The  leapfrog  procedure  was  used  to  determine  the  time  evolution  of 
the  difference  equations.  To  start  the  scheme,  a first-order  accurate 
scheme  is  used  for  the  first  time  step.  The  computer  time  required  for 
various  size  computations  "was  reported. 

At  each  time  step,  a system  of  linear  algebraic  equations,  the  finite 
difference  approximation  to  the  elliptic  equation  for  the  pressure,  must 
be  solved.  The  linear  algebraic  system  of  equations  can  be  written  in 
matrix  form,  in  which  case  the  coefficient  matrix  for  the  vector  of  unknown 
pressure  values  at  grid  points  is  large  and  sparse.  It  is  important  to 
utilize  storage  carefully  and  solve  this  algebraic  system  efficiently  and 
accurately  since  the  system  is  embedded  in  the  larger  fluid  calculation  and 
it  must  be  solved  repetitively.  A hybrid  method,  combining  the  inter- 
active algorithm  called  conjugate  gradients  with  a fast  direct  method 
for  solving  Poisson's  equation,  was  used  to  solve  the  algebraic  system. 

(it  was  found  typically  that  the  solution  of  the  algebraic  equations 
arising  on  a 31  x 31  mesh  could  be  determined  to  five  significant  figures 
after  between  two  and  five  iterations  of  the  conjugate  gradients  algorithm 
in  between  two  to  five  seconds  of  CPU  time  on  the  N.B.S.  UNIVAC  1108. ) 

Finally,  results  of  some  computations  were  presented  using  computer- 
generated plots.  A discussion  was  given  of  the  software  required  to 
generate  these  plots.  A very  brief  discussion  was  presented  of  the 
analytical  procedures  required  to  obtain  a solution  to  the  differential 
equations  describing  heating  of  a homogeneous  fluid.  Computations  of  the 
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solution  to  the  difference  equations  were  found  to  confirm  the  conclusions 
of  the  analysis.  Also  these  computations  are  found  to  yield  flow  fields 
which  behave  qualitatively  the  same  as  flow  fields  observed  in  a room 
fire.  Examples  of  computations  of  the  flow  fields  induced  by  heating  in 
a stratified  fluid  were  also  presented  and  discussed. 
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